#include "integrate1.h"
#include "common.h"
#include "common_cuda.h"

#define BLOCK_SIZE 512

__global__ void updateVel(float4* pos, float4* vel, int n, float dt)
{
	int x = blockIdx.x*BLOCK_SIZE + threadIdx.x;

	if (x < n)
	{
		float4 pos_a = pos[x];
		float3 acc = make_float3(0.0f);

		// sum total force
		int i;
		for (i = 0; i != n; ++i)
		{
			//if (i != x)
			{
				float4 pos_b = pos[i];
				acc += calcGravityAccel(&pos_a, &pos_b);
			}
		}

		// apply acceleration
		vel[x] += make_float4(dt*acc);	// 7 flops
	}
}

/*
__global__ void updatePos(float4* pos, float4* vel, int n, float dt)
{
	int x = blockIdx.x*BLOCK_SIZE + threadIdx.x;

	if (x < n)
		pos[x] += dt*vel[x];	// 3 flops
}
*/

void integrate1(float4* d_pos, float4* d_vel, int size, int iterations, float dt)
{
	dim3 dimb(BLOCK_SIZE, 1);
	dim3 dimg((size-1)/dimb.x + 1, 1);

	for (int i = 0; i != iterations; ++i)
	{
		updateVel<<<dimg, dimb>>>(d_pos, d_vel, size, dt);
		cudaThreadSynchronize();
		updatePos<<<dimg, dimb>>>(d_pos, d_vel, size, dt);
		cudaThreadSynchronize();
	}
}

void integrate1(bodyList* bl, int iterations, float dt)
{
	float4* d_pos;
	float4* d_vel;

	// allocate memory on gpu and copy data
	d_pos = (float4*) safe_cudaMalloc(bl->size*sizeof(float4));
	d_vel = (float4*) safe_cudaMalloc(bl->size*sizeof(float4));
	cudaMemcpy(d_pos, bl->pos, bl->size*sizeof(float4), cudaMemcpyHostToDevice);
	cudaMemcpy(d_vel, bl->vel, bl->size*sizeof(float4), cudaMemcpyHostToDevice);

	integrate1(d_pos, d_vel, bl->size, iterations, dt);

	// copy results back
	cudaMemcpy(bl->pos, d_pos, bl->size*sizeof(float4), cudaMemcpyDeviceToHost);
	cudaMemcpy(bl->vel, d_vel, bl->size*sizeof(float4), cudaMemcpyDeviceToHost);

	// free memory
	cudaFree(d_pos);
	cudaFree(d_vel);
}
